function exh_apeo_seq()
% 设置最大 n 的值
min_n = 17;
max_n = 17;
% 查找最大 MF 的值
search_max_mf(min_n , max_n );
% 保存结果到 CSV 文件

end
function results = search_max_mf(min_n , max_n )
% 初始化结果存储

for n = min_n:max_n
    tic;  % 记录开始时间
    % 设置 x 和 y 的范围及前缀长度
    x_k_min = max(1,floor(n / 4));
    x_k_max = floor(n / 2);
    x_prefix_length = 0;  % x 前缀的长度,由于由于循环特性

    y_k_min = max(1,floor(n / 4));
    y_k_max = ceil(n / 2)-1;
    y_prefix_length = 0;  % y 前缀的长度

    min_runs = floor(n / 2) - 4;
    max_runs = floor(n / 2) + 4;

    diff_num_best = 4 ; %保存不同解的数量

    min_peak = floor(n / 2);

    % 1. 初始化 best_params 作为结构体数组
    num_best = 4*floor(sqrt(n));  % 要保存的最佳参数数量
    
%     best_params(100) = struct('max_mf', -inf, 'best_sequences', {cell(2, 1)}, 'best_peak', 0, 'best_isl', 0, 'best_psl', 0, 'str', '');

   

%     x_all_matrix = false(0, n);  % 存储所有 x 矩阵
%     y_all_matrix = false(0, n); % 存储所有 y 矩阵
% 
%     % 生成 x 矩阵
%     parfor k = x_k_min:x_k_max
%         % 生成组合和二进制映射
%         result_map = generateCombinationsAndBinaryMap(n, k, x_prefix_length);
%         % 将二进制映射转换为矩阵
%         x_matrix = convertMapToMatrix(result_map,0,min_runs,max_runs);
%         x_all_matrix = vertcat(x_all_matrix, x_matrix); % 添加到总 x 矩阵
%     end
% 
%     % 生成 y 矩阵
%     for k = y_k_min:y_k_max
%         % 生成组合和二进制映射
%         result_map = generateCombinationsAndBinaryMap(n, k, y_prefix_length);
%         % 将二进制映射转换为矩阵
%         y_matrix = convertMapToMatrix(result_map,1,min_runs,max_runs);
%         y_all_matrix = [y_all_matrix; y_matrix];  % 添加到总 y 矩阵
%     end
    [x_all_matrix,y_all_matrix] = exh_apeo_data (n,min_runs,max_runs,x_k_min ,x_k_max , x_prefix_length,y_k_min ,y_k_max , y_prefix_length  );
    [j_max_mf_value,j_max_mf_index] = exh_apeo_gpu_fft (n,min_peak,x_all_matrix,y_all_matrix);

    [~, sortedmfs_index] = sort(j_max_mf_value(:,1), 'descend');  % 降序排序

    tmp_diff_value = 1;

    for best_index= 1:length(sortedmfs_index)
        if best_index ~= 1 && j_max_mf_value(sortedmfs_index(best_index),1) ~= j_max_mf_value(sortedmfs_index(best_index -1 ),1)
           tmp_diff_value = tmp_diff_value+1;
        end
        if best_index > num_best &&  tmp_diff_value > diff_num_best
            break;
        end
        best_params(best_index).max_mf = j_max_mf_value(sortedmfs_index(best_index,1),1);
        y=y_all_matrix(sortedmfs_index(best_index),:);
        x=x_all_matrix(j_max_mf_index(sortedmfs_index(best_index),1),:);

        R = fft_01(y).* conj(fft_01(x));  % 计算相关性
        correlation_results =  round(ifft(R, [], 2)) ;% 计算逆 FFT
        best_params(best_index).best_peak = abs(correlation_results(1));
        best_params(best_index).best_isl = sum(correlation_results(2:n).^2,2) ;
        best_params(best_index).best_psl = max(abs(correlation_results(2:n))) ;
        
        best_params(best_index).best_sequences{1} = x;
        if correlation_results(1)<0
            best_params(best_index).best_sequences{2} = ~y;
            correlation_results = ~correlation_results;
        else
            best_params(best_index).best_sequences{2} = y;
        end

         best_params(best_index).str = strcat('NO',num2str(tmp_diff_value),':( ',num2str(correlation_results(1:n)),' )') ;

    end


    elapsed_time = toc;  % 记录运行时间
    % 显示结果

    display_results(n, best_params, elapsed_time);
    save_results_to_csv(n, best_params, 'csv_apeo_seq.csv');
end
end



function [ mfs_value , mfs_index] = compute_max_mfs_limit (n,Y_FFTs_j, X_FFTs,x_all_matrix,y_seq,min_peak)
    R=zeros(size(X_FFTs,1),2*n-1);
    peaks=zeros(size(X_FFTs,1),1);
    isls=zeros(size(X_FFTs,1),1);
    mfs=zeros(size(X_FFTs,1),1);

    i = find(abs(n-2*sum(xor(y_seq , x_all_matrix),2)) > min_peak); %主峰的绝对值要大于
    
    R(i,:) = Y_FFTs_j .* conj(X_FFTs(i,:));  % 计算相关性
    correlation_results(i,:) =  round(ifft(R(i,:), [], 2));% 计算逆 FFT

    peaks(i) = correlation_results(i,1);
    isls(i) = sum(correlation_results(i,2:n).^2,2) ;
    mfs(i) = (peaks(i).^2)./isls(i)/2;

    [ mfs_value , mfs_index] = max(mfs);
end

function binary_map = generateCombinationsAndBinaryMap(n, k, prefix_length)
% 生成 n 选 k 的所有组合并返回二进制映射
combinations = nchoosek(1:n-prefix_length, k-prefix_length);  % 计算组合
binary_map = containers.Map(); % 初始化映射

% 为每个组合创建二进制字符串
num_combinations = size(combinations, 1);  % 获取组合数量
binary_string_template = zeros(1, n-prefix_length);  % 创建二进制字符串模板
prefix_array = ones(1, prefix_length);
for i = 1:num_combinations
    combination = combinations(i, :);  % 当前组合

    binary_array = binary_string_template;  % 用模板初始化二进制字符串
    binary_array(combination) = 1;  % 设置与组合对应的位置为 1
    x = [prefix_array, binary_array];  % 与前缀组合

    % 将二进制数组转换为字符串
    binary_string = num2str(x);

    % 将二进制字符串作为键存储在映射中，并给定其值
    binary_map(binary_string) = x ;
end
end

function binary_matrix = convertMapToMatrix(binary_map, bj,min_runs,max_runs)
% 将二进制映射转换为行数等于映射键计数的矩阵
keys = binary_map.keys();
binary_matrix = false(length(keys), length(binary_map(keys{1})));  % 初始化结果矩阵

row_index = 1;  % 矩阵行索引

for i = 1:length(keys)
    current_key = binary_map(keys{i});
    if ~isempty(current_key)  % 检查值是否存在

        if is_valid_runs(current_key,min_runs,max_runs),continue;end
        binary_matrix(row_index, :) = current_key;
        if bj
            binary_matrix(row_index, :) = check_equivalent_binary_sequence(current_key,binary_map); 
        end
        row_index = row_index + 1;  % 移动到下一行
    end
end

% 移除未使用的行
binary_matrix = binary_matrix(1:row_index-1, :);  % 仅保留有效行
end

function  x_fft = fft_01(array_01)
% 计算 x 和 y 的 FFT 相关性
x_seq = ones(size(array_01,1),size(array_01,2));
x_seq(array_01 == 0) = -1 ; % 将 0 转换为 -1
x_fft = fft(x_seq, size(x_seq,2)*2-1,2);  % 对 x 进行 FFT

end

function  min_seq = check_equivalent_binary_sequence(x, yset)
% 获取序列长度
n = length(x);
x = x(:)'; % 确保 x 是行向量

% 自身变换 (I)
I = x;
%if ~iskeys(num2str(I)) result=true;return;end

% 翻转 (C)
C = ~x;

% 偶数位翻转 (A)
A = x;
A(2:2:end) = ~x(2:2:end); % 翻转奇数位

% 反转 (R)
R = fliplr(x);

% 初始化等价序列集合
equivalent_sequences = {I, R, C, A}; % 包含基础变换

if mod(n, 2) == 0 % n 为偶数
    % RA: 先反转 R
    RA = R;
    RA(2:2:end) = ~RA(2:2:end); % 然后进行奇数位翻转

    % AR: 先对 A 进行奇数位翻转
    AR = A;
    AR = fliplr(AR); % 然后进行反转

    RC=R;
    RC = ~RC; % RC: C 的反转

    % CA: 先反转 C
    CA = C;
    CA(2:2:end) = ~C(2:2:end); % 在反转后的基础上进行奇数位翻转

    equivalent_sequences = [equivalent_sequences, RA, AR, RC, CA];
else % n 为奇数
    RC = ~R; % RC: R 的反转

    RA = R;
    RA(2:2:end) = ~RA(2:2:end); % RA: R 的反转

    % CA: 先反转 C
    CA = C;
    CA(2:2:end) = ~CA(2:2:end); % 在反转后的基础上进行奇数位翻转


    RCA = ~RC(2:2:end); % RCA: CA 的反转

    equivalent_sequences = [equivalent_sequences, RC, RA, CA, RCA];
end



min_seq2num = inf;
min_seq = x;


for i=1 :length(equivalent_sequences)
    seq=equivalent_sequences{i};
    keyString = num2str(seq); % 转换为字符串
    if isKey(yset, keyString)

       if bin2dec(keyString) < min_seq2num
                min_seq = seq;
                min_seq2num = bin2dec(num2str(seq));
        end
        yset(keyString) = [];
    end
end


end

function valid = is_valid_runs(seq,min_runs,max_runs)


total_runs = 1;
current_value = seq(1);

for i = 2:length(seq)
    if seq(i) ~= current_value
        total_runs = total_runs + 1;
        current_value = seq(i);
    end
end

valid = (total_runs < min_runs || total_runs > max_runs);
end
function display_results(n, params, elapsed_time)

% 显示当前 n 的结果

for i=1 :length(params)

    disp(['n = ', num2str(n)]);
    disp('最佳序列:');
    disp(params(i).best_sequences{1});  % 显示最佳序列 X
    disp(params(i).best_sequences{2});  % 显示最佳序列 Y
    disp(['主峰值: ', num2str(params(i).best_peak)]);
    disp(['MF: ', num2str(params(i).max_mf)]);
    disp(['运行时间: ', num2str(elapsed_time), '秒']);

end
disp('-----------------------');
end
function result = analyze_seq(x)
    n = length(x);
    run_seqs = [];
    current_len = 1;

    run1_max_len = 0;
    run1_count = 0;
    run0_max_len = 0;
    run0_count = 0;
    run1_len1_count = 0;

    for i = 2:n
        % 如果当前元素与前一个元素相同，增加当前长度
        if x(i) == x(i-1)
            current_len = current_len + 1;
        else
            % 将当前长度添加到序列
            run_seqs = [run_seqs current_len];
            if x(i-1)==1 && current_len==1
                run1_len1_count = run1_len1_count+1;
            end

            % 处理当前序列的计数和最大长度
            if x(i-1) == 1
                if current_len > run1_max_len
                    run1_max_len = current_len;
                    run1_count = 1;
                elseif current_len == run1_max_len
                    run1_count = run1_count + 1;
                end
            else
                if current_len > run0_max_len
                    run0_max_len = current_len;
                    run0_count = 1;
                elseif current_len == run0_max_len
                    run0_count = run0_count + 1;
                end
            end

            % 重置当前长度
            current_len = 1;
        end
    end

    % 处理最后一个序列
    run_seqs = [run_seqs current_len];

    if x(end) == 1
        if current_len > run1_max_len
            run1_max_len = current_len;
            run1_count = 1;
        elseif current_len == run1_max_len
            run1_count = run1_count + 1;
        end
    else
        if current_len > run0_max_len
            run0_max_len = current_len;
            run0_count = 1;
        elseif current_len == run0_max_len
            run0_count = run0_count + 1;
        end
    end

    % 过滤运行序列
    run_seqs = run_seqs(run_seqs > 0);

    % 计算差异
    run_count = length(run_seqs);
    
    run_count_diff = floor(n / 2) - run_count;
    run1_max_len_diff = run1_max_len - floor(sqrt(n));
    run0_max_len_diff = run0_max_len - floor(sqrt(n));

    % 创建结果结构
%    result.run_seqs = run_seqs;
    result.k = sum(x);

    result.k_diff = floor(n / 2) - result.k;

    result.run_count = run_count;
    result.run_count_diff = run_count_diff;
    result.run1_count = run1_count;
    result.run1_count_diff = floor(n / 4) - run1_count;
    result.run0_count = run0_count;
    result.run0_count_diff = floor(n / 4) - run0_count;
    result.run1_max_len = run1_max_len;
    result.run1_max_len_diff = run1_max_len_diff;
    result.run0_max_len = run0_max_len;
    result.run0_max_len_diff = run0_max_len_diff;
    result.run1_len1_count = run1_len1_count;
    result.run1_len1_count_diff = floor(result.k / 2) - run1_len1_count;
end
